
// c++ includes
#include <tr1/functional>
#include <ostream>

#include "libmesh.h"

#include "neutronics_system.h"
#include "wrapper.h"

#include "coupling_matrix.h"
#include "string_to_enum.h"

#include "mesh.h"
#include "dof_map.h"
#include "explicit_system.h"
#include "slepc_eigen_solver.h"

#include "utility.h"
//#include "petsc_vec_utils.h"
#include "string_to_enum.h"

// Define the Finite Element object.
#include "fe.h"

// Define Gauss quadrature rules.
#include "quadrature_gauss.h"

#include "problem.h"


PerfLog NeutronicsSystem::logger("EigenSolver") ;

#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::constructor"
NeutronicsSystem::NeutronicsSystem (libMesh::MeshBase& o_mesh, const Problem& pb, 
									const unsigned int e_groups,
									const std::string s_order,
									const std::string s_family) :
		EquationSystems(o_mesh),
		_system (EquationSystems::add_system<ExplicitSystem>("NeutronicsEigenProblem")),
		_egroups (e_groups), _problem(&pb), 
		_recompute_group_blocks (true),
		_is_eigen_problem (true),
		_is_eigen_problem_mfree (false),
		_only_diagonal_loss_matrix (false),
		_is_gain_matrix_created (false),
		_is_loss_matrix_created (false),
		_only_block_diagonal_loss_matrix (false),
		_keff (1.0),
		_eigen_shift (0.0)
{
	Order _order = Utility::string_to_enum<Order>(s_order) ;
	std::stringstream out;

	for (unsigned int ieg=0; ieg < e_groups; ieg++)
	{
		std::string flux(ENERGY_FLUXES) ;
		out << ieg ;
		flux.replace(flux.find("*"), 1, out.str()) ;

		// loop over e_groups and add variables to _system
		_system.add_variable( flux, Utility::string_to_enum<Order>(s_order), Utility::string_to_enum<FEFamily>(s_family) );
		out.str("") ;
	}

	// set the default quadrature order
	//_default_quad_order = static_cast<Order>(2*(int)_order + 1) ;

	_initialized = false ;
}

#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::destructor"
NeutronicsSystem::~NeutronicsSystem ()
{
	// Clear data
	this->clear();
}

#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::initialize_eigen_context"
void NeutronicsSystem::initialize_eigen_context(bool is_mfree, EigenProblemType eigen_type)
{
	_is_eigen_problem = true ;

	_is_eigen_problem_mfree = is_mfree ;

	_eigen_problem_type = eigen_type ;

	// create the linear solver for Loss operator
	PetscErrorCode ierr = 0 ;
	// Set the block KSP options
	ierr = KSPCreate(PETSC_COMM_WORLD, &_loss_matrix_solver) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

	// we will take the options if it is prefixed with "lossoptr_"
	ierr = KSPAppendOptionsPrefix(_loss_matrix_solver, "lossoptr_") ;

	// since the bigger block matrices is unsymmetric, let us use BICGSTAB as default solver
	ierr = KSPSetType(_loss_matrix_solver, KSPGMRES) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

	// let us set LU as the preconditioner only if specified in command line
	PC lossmatpc ;
	ierr = KSPGetPC(_loss_matrix_solver, &lossmatpc);CHKERRABORT(PETSC_COMM_WORLD,ierr);
	ierr = PCSetType(lossmatpc, PCSHELL) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
	ierr = PCShellSetApply(lossmatpc, Wrapper::precondition_loss_matrix_solve) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
	ierr = PCShellSetContext(lossmatpc, this);CHKERRABORT(PETSC_COMM_WORLD,ierr);
	ierr = PCShellSetName(lossmatpc,"block-gauss-seidel");CHKERRABORT(PETSC_COMM_WORLD,ierr);

	ierr = KSPSetFromOptions(_loss_matrix_solver) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

}


//#undef __FUNCT__
//#define __FUNCT__ "NeutronicsSystem::create_loss_matrix"
//void NeutronicsSystem::create_loss_matrix() 
//{
//	if (!_is_eigen_problem_mfree && !_is_loss_matrix_created)
//	{
//		_is_loss_matrix_created = true ;
//
//		this->assemble_loss_matrix(/*_loss_mat*/) ;
//	}
//
//	return _loss_mat ;
//}
//
//#undef __FUNCT__
//#define __FUNCT__ "NeutronicsSystem::create_gain_matrix"
//void NeutronicsSystem::create_gain_matrix() 
//{
//	if (!_is_eigen_problem_mfree && !_is_gain_matrix_created)
//	{
//		_is_gain_matrix_created = true ;
//
//		this->assemble_gain_matrix(/*_gain_mat*/) ;
//	}
//}


#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::get_loss_matrix"
PetscMatrix<Real>* NeutronicsSystem::get_loss_matrix() const
{
	return _loss_mat ;
}

#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::get_gain_matrix"
PetscMatrix<Real>* NeutronicsSystem::get_gain_matrix() const
{
	return _gain_mat ;
}


#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::solve_eigen_problem"
Real NeutronicsSystem::solve_eigen_problem(Real tol, Real maxit)
{
    PetscErrorCode ierr = 0 ;

    SlepcEigenSolver<Real> eig_solver ;
    Real relerr = 1.0 ;

	_system.update() ;

    // recompute the gain and loss matrices and compute fundamental eigenmode.
	ierr = this->assemble_loss_matrix () ;CHKERRQ(ierr);
	ierr = this->assemble_gain_matrix () ;CHKERRQ(ierr);
	
    eig_solver.set_eigensolver_type(KRYLOVSCHUR) ;
    eig_solver.set_eigenproblem_type(GNHEP) ;
    eig_solver.set_position_of_spectrum(LARGEST_MAGNITUDE) ;
	
	std::cout << "\nabout to solve the eigenproblem." ;
	
    std::pair<unsigned int, unsigned int> eig_result = eig_solver.solve_generalized(*_gain_mat,
                       *_loss_mat,
                       1,	// we want only fundamental mode
                       30,	// set default subspace size
                       tol,	// can make this user specified later
                       maxit);	// for Krylov-Schur, this should be plenty.

	std::cout << "\nsolved eigenproblem. now normalizing" ;

    if (eig_result.first >= 1)	// we asked for only fundamental mode; see if we converged
    {
		// make sure we create a vector which contains only the flux unknowns
		// TODO: eliminate precursor equation through Gaussian elimination so that
		// only flux values are in the primary system
        PetscVector<Real>& flx_soln = this->get_solution() ;
		
        std::pair<Real,Real> eig = eig_solver.get_eigenpair(0, flx_soln) ;	flx_soln.close() ;
        _keff = eig.first ;
		relerr = eig_solver.get_relative_error(0) ;
		PetscPrintf(PETSC_COMM_WORLD,"\nNeutronics eigenproblem: Keff = %14.15G, Relative error = %G, Iterations = %D\n", _keff, relerr, eig_result.second);
		
		Parallel::barrier() ;
		
        bool normalize_solution = true ;
        // normalzie fundamental mode
        if (normalize_solution)
        {
            PetscScalar maxv = 0., minv = 0. ;	int ind = -1 ;

			// a bad hack -> to offset negative values  at strong material interfaces
			// this causes problem during a transient
			ierr = VecAbs(flx_soln.vec()) ;CHKERRQ(ierr);
			
            ierr = VecNorm(flx_soln.vec(), NORM_INFINITY, &maxv) ;CHKERRQ(ierr);
            ierr = VecScale(flx_soln.vec(), 1.0/maxv) ;CHKERRQ(ierr);	flx_soln.close() ;
            ierr = VecMax(flx_soln.vec(), &ind, &maxv) ;CHKERRQ(ierr);

            if (fabs(1.0-maxv) > 1e-10)	// possibly negative maxima
            {
                ierr = VecScale(flx_soln.vec(), -1.0) ;CHKERRQ(ierr);	flx_soln.close() ;
            }

            // ierr = VecMax(flx_soln.vec(), &ind, &maxv) ;CHKERRQ(ierr);
            // ierr = VecMin(flx_soln.vec(), &ind, &minv) ;CHKERRQ(ierr);
            // PetscPrintf(PETSC_COMM_WORLD, "\nMaximum value = %G and Minimum value = %G\n", maxv, minv) ;
        }

		this->update() ;

		/// Perform our own check to compute the residual of the eigenvalue problem
		// Vec pvtemp, pvres ;
		// ierr = VecDuplicate(flx_soln.vec(), &pvtemp) ;CHKERRQ(ierr);
		// ierr = VecDuplicate(flx_soln.vec(), &pvres) ;CHKERRQ(ierr);
        // PetscVector<Real> temp(pvtemp) ;
        // PetscVector<Real> res(pvres) ;
		// ierr = MatMult(lossmat->mat(), flx_soln.vec(), res.vec()) ;	res.close() ;
		// ierr = MatMult(gainmat->mat(), flx_soln.vec(), temp.vec()) ;	temp.close() ;
		// ierr  = VecAXPY(res.vec(), -1.0/_keff, temp.vec()) ;			res.close() ;
		// relerr = res.l2_norm() ;
		// temp.clear() ; res.clear() ;
		// ierr = VecDestroy(&pvtemp) ;CHKERRQ(ierr);
		// ierr = VecDestroy(&pvres) ;CHKERRQ(ierr);
		
		// soln.insert(flx_soln, flx_rows) ;	soln.close() ;
		
		// commenting out the compute_power call 
		//params.set<bool>("update_power_solution") = false ;
        //Real total_power = this->compute_power(this->get_power(), params) ;

		// PetscPrintf(PETSC_COMM_WORLD,"\nTime=%G, Total_power=%5.12E", _ctx->time, total_power) ;
		
		// this->normalize_flux_values() ;
	
    }
    else
    {
        std::cout << "\n  Baaah .. Not sure if I converged. Pausing. " ;
		std::cout << "\n  Run again with -eps_monitor to see if converged. " ;
        std::cin.get() ;
    }

    return relerr ;
}

#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::clear"
void NeutronicsSystem::clear ()
{
	if (_initialized)
	{
		PetscErrorCode ierr = 0;

		// we should probably clear the memory allocated for the eigenvalue calculation.
		// check if eigencontext = true and then deallocate matrices, unless there is a 
		// use for it as a preconditioner !
		if (_is_eigen_problem_mfree)
		{
			Mat A = _loss_mat->mat() ;
			delete _loss_mat ;
			ierr = MatDestroy(&A) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
		
			Mat B = _gain_mat->mat() ;
			delete _gain_mat ;
			ierr = MatDestroy(&B) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
		}
		// Destroy the loss matrix solver.
		ierr = KSPDestroy(&_loss_matrix_solver) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
		
		_is_eigen_problem = false ;
		_is_eigen_problem_mfree = false ;
		_only_diagonal_loss_matrix = false ;
		_is_gain_matrix_created = false ;    // can reuse this matrix for preconditioning
		_is_loss_matrix_created = false ;    // can reuse this matrix for preconditioning
		_only_block_diagonal_loss_matrix = false ;
	
		if(_group_block_vectors != NULL)
		{
			// delete our group block matrices and solvers used for preconditioning
			for (unsigned int ie=0; ie < _egroups; ie++)
			{
				delete _group_block_vectors[ie] ;
				_group_block_vectors[ie] = NULL ;

				for (unsigned int je=0; je <= ie; je++)
				{
					delete _group_block_matrices[ie*_egroups+je] ;
					_group_block_matrices[ie*_egroups+je] = NULL ;
				}

				//delete _group_block_solvers[ie] ;
				ierr = KSPDestroy(&_group_block_solvers[ie]) ;	CHKERRABORT(PETSC_COMM_WORLD,ierr);
			}
		}
		
		// clear all cached variables
		flx_var.clear() ;

		_initialized = false ;
	}
}

#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::init"
void NeutronicsSystem::init ()
{
	if (!_initialized)
	{
		// if we have delayed groups and we are doing the eigenvalue problem
		// do not allocate the matrix yet.
		//_system.delay_matrix_initialization() = (_is_eigen_problem && (this->precursor_groups() > 0)) ;
		
		// initialize parent data
		EquationSystems::init();

		std::stringstream out;
		// compute all the variable numbers
		{
			flx_var.resize(_egroups) ;

			for (unsigned int ieg=0; ieg < _egroups; ieg++)
			{
				std::string flux(ENERGY_FLUXES) ;
				out << ieg ;
				flux.replace(flux.find("*"), 1, out.str()) ;

				// loop over e_groups and find variable number
				flx_var[ieg] = (unsigned int)_system.variable_number( flux );
				
				out.str("") ;
			}

		}

		// variable numbers computed

		PetscErrorCode ierr = 0 ;
	
		// Allocate a single block/group which can be used in Block-Jacobi and/or Gauss-seidel preconditioning
		unsigned int nvars = _system.n_vars() ;

		const DofMap& dof_map = _system.get_dof_map() ;
		
		const std::vector<unsigned int>& n_nz = dof_map.get_n_nz();
		const std::vector<unsigned int>& n_oz = dof_map.get_n_oz();

		// // // Ripped off the following part from libMesh's DofMap code
		// CouplingMatrix dof_coupling (nvars) ;
		
		// for (unsigned int ig=0; ig < _egroups; ig++)
			// for (unsigned int jg=0; jg < _egroups; jg++)
				// dof_coupling(ig,jg) = 1 ;
		
		// bool implicit_neighbor_dofs = false, need_full_sparsity_pattern = false ;
		// // We can compute the sparsity pattern in parallel on multiple
		// // threads.  The goal is for each thread to compute the full sparsity
		// // pattern for a subset of elements.  These sparsity patterns can
		// // be efficiently merged in the SparsityPattern::Build::join()
		// // method, especially if there is not too much overlap between them.
		// // Even better, if the full sparsity pattern is not needed then
		// // the number of nonzeros per row can be estimated from the
		// // sparsity patterns created on each thread.
		// SparsityPattern::Build sp (this->get_mesh(),
					 // dof_map,
					 // &dof_coupling,
					 // implicit_neighbor_dofs,
					 // need_full_sparsity_pattern);

		// Threads::parallel_reduce (ConstElemRange (this->get_mesh().active_elements_begin(),
							// this->get_mesh().active_elements_end()), sp);

		// // steal the n_nz and n_oz arrays from sp and use it to build our block
		// // matrix because this is more efficient.
		// // (sp.n_nz), (sp.n_oz)
		
		// const std::vector<unsigned int>& n_nz = sp.n_nz;
		// const std::vector<unsigned int>& n_oz = sp.n_oz;
		
		if (_is_eigen_problem )
		{
			PetscErrorCode ierr = 0;

			if(!_is_eigen_problem_mfree)
			{
				// Allocate a single block/group which can be used in Block-Jacobi and/or Gauss-seidel preconditioning

				unsigned int nglobdofs = _system.n_dofs() ;
				unsigned int nlocdofs = _system.n_local_dofs() ;
				
				std::cout << "\nNvars = " << nvars << "\tEnergy groups = " << _egroups ;
				std::cout << "\nGlobal: ndofs = " << _system.n_dofs() << "\tlocal_ndofs = " << _system.n_local_dofs() ;

				// Create and allocate correctly for both Gain and Loss matrices
				_gain_mat = new PetscMatrix<Real>() ;
				_gain_mat->attach_dof_map(dof_map) ;
				_gain_mat->init() ;

				_loss_mat = new PetscMatrix<Real>() ;
				_loss_mat->attach_dof_map(dof_map) ;
				_loss_mat->init() ;
				
				// std::cout << "Eigen: ndofs = " << nglobdofs << "\tlocal_ndofs = " << _loss_mat->row_stop()-_loss_mat->row_start() << ", " << _gain_mat->row_stop()-_gain_mat->row_start() ;
				
				// _loss_mat->init(nglobdofs, nglobdofs, nlocdofs, nlocdofs, n_nz[0], n_oz[0]) ;
			}
			else
			{
				// NOTE: not sure if n_dofs is available yet; we might have to post-pone this operation
				// till during init(). need check this further for mat-free runs..
				int numdofs = _system.n_dofs() ;
				int locdofs = _system.n_local_dofs() ;

				Mat shellGain, shellLoss ;

				ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)this, &shellLoss);CHKERRABORT(PETSC_COMM_WORLD,ierr);
				ierr = MatSetFromOptions(shellLoss) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
				ierr = MatShellSetOperation(shellLoss, MATOP_MULT,(void(*)())Wrapper::multiply_loss_matrix);CHKERRABORT(PETSC_COMM_WORLD,ierr);

				ierr = MatCreateShell(PETSC_COMM_WORLD,locdofs,locdofs,numdofs,numdofs,(void*)this, &shellGain);CHKERRABORT(PETSC_COMM_WORLD,ierr);
				ierr = MatSetFromOptions(shellGain) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);
				ierr = MatShellSetOperation(shellGain, MATOP_MULT,(void(*)())Wrapper::multiply_gain_matrix);CHKERRABORT(PETSC_COMM_WORLD,ierr);

				_gain_mat = new PetscMatrix<Real> (shellGain) ;
				_loss_mat = new PetscMatrix<Real> (shellLoss) ;
			}
			
			// we might need this for eigenvalue problem solution using Newton
			// as a preconditioner; need a better flag for this purpose
			// if (_only_block_diagonal_loss_matrix)
			{
				// NOTE: This should also be called even when loss_mat is shell.
				ierr = KSPSetOperators(_loss_matrix_solver, _loss_mat->mat(), _loss_mat->mat(), SAME_PRECONDITIONER) ;
				CHKERRABORT(PETSC_COMM_WORLD,ierr);
				
				const unsigned int nglobdofs = _loss_mat->m() ;
				const unsigned int nlocdofs = _loss_mat->row_stop()-_loss_mat->row_start() ;
					
				// allocate for our block matrices and vectors
				_group_block_matrices = new PetscMatrix<Real>*[_egroups*_egroups] ;
				_group_block_vectors = new PetscVector<Real>*[_egroups] ;
				_group_block_solvers = new KSP[_egroups] ;

				// create the block vectors and matrices according to dof-map sparsity
				for (unsigned int ie=0; ie < _egroups; ie++)
				{
					_group_block_vectors[ie] = new PetscVector<Real>(nglobdofs, nlocdofs) ;

					// NOTE: the preconditioner is only block-gauss-seidel. so we only need diagonal and
					// lower triangular blocks to exist. for future preconditioners, alternate variations
					// can be created.
					for (unsigned int je=0; je <= ie; je++)
					{
						_group_block_matrices[ie*_egroups+je] = new PetscMatrix<Real>() ;
						_group_block_matrices[ie*_egroups+je]->init(nglobdofs, nglobdofs, nlocdofs, nlocdofs, n_nz[0], n_oz[0]) ;
					}
				}

				// create KSP solvers for all groups and associate with block-diagonal group matrices
				for (unsigned int ie=0; ie < _egroups; ie++)
				{
					// Set the block KSP options
					ierr = KSPCreate(PETSC_COMM_WORLD, &_group_block_solvers[ie]) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

					std::string flux(ENERGY_FLUXES) ;
					out << "group_" << ie << "_" ;

					// we will take the options if it is prefixed with "group_ie_"
					ierr = KSPAppendOptionsPrefix(_group_block_solvers[ie], out.str().c_str()) ;

					// set the matrix references to KSP object
					ierr = KSPSetOperators(_group_block_solvers[ie], _group_block_matrices[ie*_egroups+ie]->mat(), _group_block_matrices[ie*_egroups+ie]->mat(), SAME_NONZERO_PATTERN) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

					ierr = KSPSetTolerances(_group_block_solvers[ie], 1e-10, 1e-14, PETSC_DEFAULT, 500) ; CHKERRABORT(PETSC_COMM_WORLD,ierr);

					ierr = KSPSetInitialGuessKnoll(_group_block_solvers[ie], PETSC_TRUE) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

					// since the block matrices are symmetric, let us use PCG as default solver
					ierr = KSPSetType(_group_block_solvers[ie], KSPCG) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

					// let us set ICC as the default preconditioner
					PC grpmatpc ;
					ierr = KSPGetPC(_group_block_solvers[ie], &grpmatpc);CHKERRABORT(PETSC_COMM_WORLD,ierr);
					ierr = PCSetType(grpmatpc, PCICC) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

					ierr = KSPSetFromOptions(_group_block_solvers[ie]) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

					out.str("") ;
				}
			}
		}
		else
		{
			_group_block_vectors = NULL ;
			_group_block_matrices = NULL ;
			_group_block_solvers = NULL ;
		}

		_initialized = true ;
	}
}

#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::reinit"
void NeutronicsSystem::reinit()
{
	// first call our base class reinit() member
	EquationSystems::reinit() ;

	if (_is_eigen_problem )
	{
		PetscErrorCode ierr ;
		// first destroy all the group vectors
		{
			// delete our group block matrices and solvers used for preconditioning
			for (unsigned int ie=0; ie < _egroups; ie++)
			{
				delete _group_block_vectors[ie] ;

				for (unsigned int je=0; je <= ie; je++)
				{
					delete _group_block_matrices[ie*_egroups+je] ;
				}

				//delete _group_block_solvers[ie] ;
				ierr = KSPDestroy(&_group_block_solvers[ie]) ;	CHKERRABORT(PETSC_COMM_WORLD,ierr);
			}
		}
		
		// now reallocate all the group vectors
		{
			// Allocate a single block/group which can be used in Block-Jacobi and/or Gauss-seidel preconditioning
			unsigned int nvars = _system.n_vars() ;
			unsigned int nlocdofs = _system.n_local_dofs()/nvars ;
			unsigned int nglobdofs = _system.n_dofs()/nvars ;

			_group_block_matrices = new PetscMatrix<Real>*[_egroups*_egroups] ;
			_group_block_vectors = new PetscVector<Real>*[_egroups] ;
			_group_block_solvers = new KSP[_egroups] ;

			const std::vector<unsigned int>& n_nz = _system.get_dof_map().get_n_nz();

			// create the block vectors and matrices according to dof-map sparsity
			for (unsigned int ie=0; ie < _egroups; ie++)
			{
				_group_block_vectors[ie] = new PetscVector<Real>(nglobdofs, nlocdofs) ;

				// NOTE: the preconditioner is only block-gauss-seidel. so we only need diagonal and
				// lower triangular blocks to exist. for future preconditioners, alternate variations
				// can be created.
				for (unsigned int je=0; je <= ie; je++)
				{
					_group_block_matrices[ie*_egroups+je] = new PetscMatrix<Real>() ;
					_group_block_matrices[ie*_egroups+je]->init(nglobdofs, nglobdofs, nlocdofs, nlocdofs, n_nz[0], n_nz[0]) ;
				}
			}

			std::stringstream out;
			// create KSP solvers for all groups and associate with block-diagonal group matrices
			for (unsigned int ie=0; ie < _egroups; ie++)
			{
				// Set the block KSP options
				ierr = KSPCreate(PETSC_COMM_WORLD, &_group_block_solvers[ie]) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

				std::string flux(ENERGY_FLUXES) ;
				out << "group_" << ie << "_" ;

				// we will take the options if it is prefixed with "group_ie_"
				ierr = KSPAppendOptionsPrefix(_group_block_solvers[ie], out.str().c_str()) ;

				// set the matrix references to KSP object
				ierr = KSPSetOperators(_group_block_solvers[ie], _group_block_matrices[ie*_egroups+ie]->mat(), _group_block_matrices[ie*_egroups+ie]->mat(), SAME_NONZERO_PATTERN) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

				ierr = KSPSetTolerances(_group_block_solvers[ie], 1e-10, 1e-14, PETSC_DEFAULT, 500) ; CHKERRABORT(PETSC_COMM_WORLD,ierr);

				ierr = KSPSetInitialGuessKnoll(_group_block_solvers[ie], PETSC_TRUE) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

				// since the block matrices are symmetric, let us use PCG as default solver
				ierr = KSPSetType(_group_block_solvers[ie], KSPCG) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

				// let us set ICC as the default preconditioner
				PC grpmatpc ;
				ierr = KSPGetPC(_group_block_solvers[ie], &grpmatpc);CHKERRABORT(PETSC_COMM_WORLD,ierr);
				ierr = PCSetType(grpmatpc, PCICC) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

				ierr = KSPSetFromOptions(_group_block_solvers[ie]) ;CHKERRABORT(PETSC_COMM_WORLD,ierr);

				out.str("") ;
			}
		}
		
	}
	
}



#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::__neutronics_exact_value"
// Exact solution function prototype.
Real __neutronics_exact_value (const Point& p,
        		      const Parameters& params,   // EquationSystem parameters, not needed
        		      const std::string&,  // sys_name, not needed
        		      const std::string& /*var_name*/) // unknown_name, not needed);
{
	const Problem* problem = params.get<const Problem*>("problem") ;

//	int index = var_name.find("flux_") ;
//	std::string varname(var_name);
//	std::istringstream i ( varname.replace(index, 5, "") ) ;
//	int egroup_val = 0;
//	// if (! (i >> egroup_val))
//		// libmesh_error() ;
//	varpars.set<unsigned int>("e_group") = static_cast<unsigned int>(egroup_val) ;

	return problem->ExactSolution(p/*, var_name*/) ;
}


#undef __FUNCT__
#define __FUNCT__ "NeutronicsSystem::form_initial_guess"
PetscErrorCode NeutronicsSystem::form_initial_guess(PetscVector<Real>& xsoln)
{
	const MeshBase& mesh = _system.get_mesh();
	const int dim = mesh.mesh_dimension();

	std::cout << "\nNumber of variables in primary system = " << _system.n_vars() ;

	this->parameters.set<const Problem*>("problem") = _problem ;
	// give an empty parameter object to make LibMesh happy;
	// _system.project_solution(__neutronics_exact_value, NULL, this->parameters);
	_system.project_vector(__neutronics_exact_value, NULL, this->parameters, xsoln);

	return 0 ;
}
